Functional renormalization group for nonequilibrium quantum many-body problems 
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We extend the concept of the functional renormahzation for quantum many-body problems to 
nonequilibrium situations. Using a suitable generating functional based on the Keldysh approach, we 
derive a system of coupled differential equations for the m-particle vertex functions. The approach 
is completely general and allows calculations for both stationary and time-dependent situations. 
As a specific example we study the stationary state transport through a quantum dot with local 
Coulomb correlations at finite bias voltage employing two different truncation schemes for the infinite 
hierarchy of equations arising in the functional renormalization group scheme. 

PACS numbers: 71.27.-|-a,73.21.La,73.23.-b 



I. INTRODUCTION 

The reliable calculation of physical properties of 
interacting quantum mechanical systems presents a 
formidable task. Typically, one has to cope with the in- 
terplay of different energy-scales possibly covering several 
orders of magnitude even for simple situations. Approxi- 
mate tools like perturbation theory, but even numerically 
exact techniques can usually handle only a restricted win- 
dow of energy scales and are furthermore limited in their 
applicability by the approximations involved or the com- 
putational resources available. In addition, due to the 
divergence of certain classes of Feynman diagrams, some 
of the interesting many-particle problems cannot be tack- 
led by straight forward perturbation theory. 

The situation becomes even more involved if one is in- 
terested in properties off equilibrium, in particular time- 
dependent situations. A standard approach for such 
cases is based on the Keldysh formalism^ for the time 
evolution of Green functions, resulting in a matrix struc- 
ture of propagators and self-energies. This structure is a 
direct consequence of the fact that in nonequilibrium we 
have to calculate averages of operators taken not with re- 
spect to the ground state but with respect to an arbitrary 
state. Therefore, the Gell-Mann and Low theorem- is not 
valid any more. Other approaches attempt to treat the 
time evolution of the nonequilibrium system numerically, 
for example using the density matrix renormalization 
groupS, the numerical renormalization group (NRG)^'^ 
or the Hamiltonian based flow-equation method.® 

The Keldysh technique shows a big flexibility and 
one therefore can find a wide range of applications 
such as transport through atomic, molecular and nano 
devices under various conditions,'' systems of atoms 
interacting with a radiation field in contact with a 
bathf^i^ or electron-electron interaction in a weakly ion- 
ized plasmaJ^ 

One powerful concept to study interacting many- 
particle systems is the rather general idea of the renor- 
malization group^^ (R-G), which has also been applied 
to time-dependent and stationary nonequilibrium situa- 
tions recently ^^i^^'^'^i^^'^® In the RG approach one usu- 



ally starts from high energy scales, leaving out possible 
infrared divergences and works ones way down to the 
desired low-energy region in a systematic way. However, 
the precise definition of "systematic way" does in general 
depend on the problem studied. 

In order to resolve this ambiguity for interacting 
quantum mechanical many-particle systems in equilib- 
rium, two different schemes attempting a unique, prob- 
lem independent prescription have emerged during the 
past decade. One is Wegner's Hamiltonian based fiow- 
equation technique fiii^ the second a field theoretical ap- 
proach, which we want to focus on in the following. This 
approach is based on a functional representation of the 
partition function of the system and has become known 
as functional renormalization group (fRG)ii^i22i2iiS^. 

A detailed description of the various possible imple- 
mentations of the fRG and its previous applications in 
equilibrium can be found e.g. in Refs. [23| and [2^. In the 
present work we extend the fRG to nonequilibrium by for- 
mulating the problem on the real instead of the imaginary 
time axis using the functional integral representation of 
the action on the Keldysh contour. Within a diagram- 
matic approach a similar set of fRG flow equations has 
already been derived by Jakobs and Schoeller and applied 
to study nonlinear transport through one-dimensional 
correlated electron systemsi ^^'^® We believe that this 
method will enable us to treat a variety of nonequilib- 
rium problems within a scheme which is well established 
in equilibrium and in contrast to other approaches is com- 
paratively modest with respect to the computer resources 
required. Our framework for nonequilibrium will turn 
out to be sufficiently general to allow for a treatment 
of systems disturbed by arbitrary external flelds (bias 
voltage, laser field, etc.), which can be constant or time- 
dependent. 

For classical many-body problems the fRG (also 
called nonperturbative renormalization group) was used 
to study nonequilibrium transitions between stationary 
states*22 

As a simple but nontrivial application to test the po- 
tential and weakness of our implementation of a nonequi- 
librium fRG we choose the single impurity Anderson 
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FIG. 1: (color online) Keldysh contour 

model (SIAM).^® This model represents the paradigm for 
correlation effects in condensed matter physics and is at 
the heart of a large range of experimenta l^^'^^i'^^'^^i'^^'^'^i'^^ 
and theoretical investigations i^^^^Z, It furthermore is the 
standard model for the description of the transport prop- 
erties of interacting single-level quantum dots. 

The paper is organized as follows. In Sec. [IT] and the 
appendices we extend the equilibrium fRG to treat gen- 
eral nonequilibrium quantum many-body problems. As 
an example we study the finite bias, stationary transport 
through an Anderson impurity in Sec. IIIII Within the 
nonequilibrium extension we obtain a system of coupled 
tensor equations, which represent the flow of the differ- 
ent components (in the Keldysh space) of the self-energy 
and the vertex function. We discuss two approximations 
derived from the fRG scheme that have been successfully 
applied to the equilibrium situation. '^^'^^ A summary and 
outlook in Sec. IIVI concludes the paper. 



II. FUNCTIONAL RG IN NONEQUILIBRIUM 

We start from the standard definition of the two-time 
Green function 



(1) 



in the interaction picture, where ^ comprises a set of 
single-particle quantum numbers and time t. The time 
evolution operators are given as 



S" = Texp i j Vi{t)dt 



(2) 



where t, t' € Ck_ , the anti time-ordered Green function 

= -z9{t t'mai^Ho) - mt' - mnme)), (s) 

with t, t' e Ck^ and 

= -imOi'HO) :t' (^CK^,t(^CK_ (6) 

where C = -1-1 for bosons and C = — 1 for fermions. Con- 
cerning the Keldysh indices we here follow the notation 
of Ref.liO. G—{^',0 and G++(C',0 take into account 
the excitation spectrum while G"' (C',?) and G ^ 
describe the thermodynamic state of the system. Only 
three of the Green functions are independent and one 
commonly introduces the linear combinations 



G«(C',0 - 0{t'-t)[G' 
G^(^',0 = 0it-t')[G- 
G^(^',0 = G-+(e',e)4 



-ie,0 + CG-+{e,0] (8) 
+ (e',0 + CG+-(^',0] (9) 
G+-(e',0 , (10) 



named retarded, advanced and Keldysh component, re- 
spectively. 

To derive a nonequilibrium fRG scheme for an inter- 
acting many-body problem along the lines of Refs. and 
[2^ one needs a formulation that allows to express the m- 
particle vertices as functional derivatives of a generating 
functional. Here we use the approach by Kamcnccv.'*- 
To set up a functional integral representation of the gen- 
erating functional we define the matrix 



G 



G- 
G^ 



G- 
G^ 



and the short hand notation 

00 



Texp ■ 



Vi{t)dt 



(3) 



with T the usual time ordering and T the anti time- 
ordering operator. The interaction term Vi(t) is arbi- 
trary, including possible explicit time dependence. 

In a nonequilibrium situation, the propagation of the 
system from — cx) ^ 00 is not any more equivalent to 
the propagation from 00 —f —00, i.e. one has to distin- 
guish whether the time arguments in Eq. ([T|) belong to 
the former or latter4Si^i^ This scheme is usually de- 
picted by the Keldysh double-time contour shown in Fig. 
[1] where Gk_ represents the propagation —00 — > 00 (up- 
per branch of the Keldysh contour) and Ck+ the propa- 
gation 00 — *■ —00 (lower branch of the Keldysh contour). 
Consequently, one has to introduce four distinct propa- 
gators, namely the time-ordered Green function 

G-{^',o = -w{eHHm 

= -i9{t' - tma^Ho) - mt - i')(^no^(r)>, (4) 



where 



is a spinor of fields (Grassmann for fermions or complex 
for bosons) with i/)- having a time from the upper branch 
of the Keldysh contour and ■0+ a time from the lower. 
Later it will also prove useful to Fourier transform from 
time t to frequency uj. One then has to replace i in ^ 
by oj. The integrals over ^ and ^' stand for summations 
over the quantum numbers and integrations over time or 
frequency. The following steps can be performed with ^ 
either containing time or frequency. The generalization 
of the functional integral representation of the partition 
function to nonequilibrium is^ 
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Go 



^ - J^int im, W) \ , So = / PV^Vexp <^ ^, 



Go 



(11) 



The matrix Go denotes the propagator of the noninteracting part of the many-body problem and 5*1111 is the (arbitrary) 
interaction term. 

To construct a generating functional for m-particle Green functions, one introduces external source fields 77 and fj 



The (connected) m-particle Green function can then be obtained by taking functional derivatives 



Gr 



V' -*5int(W,M)- (V^,^) -(^,^) 



(12) 



?7— 77— 



with 



wMM,W) = in[>v(M,M)] 



(13) 



(14) 



. . ^ . _ . . . (c) 

The derivatives in Eq. (I13p are taken with respect to the spinors and 77, i.e. the resulting Green function Gm is a 
tensor of rank 2m in the Keldysh indices. In the following we adopt the notation that quantities with explicit index 
m are tensors of rank 2m; without index m we denote the propagator and self-energy, both carrying a hat to point 
out their matrix structure. 
Introducing the fields 



we can perform a Legendre transformation 

r im, m) = -w^ m, m) - (0, v) - (^, 0) + ( ^, 



Go 



(15) 



(16) 



to the generating functional of the one-particle irreducible vertex functions 7™ (tensors of rank 2m in the Keldysh 
indices) 



-r(W,m) 



(17) 



0=0=0 



In contrast to the usual definition'*'* of T for convenience 

-1 



(see Appendix A) we added a term 



Go 



in Eq. 



(|16p . The relations between the Gm and 7™ m imaginary 
time can be found in text books^ and can straightfor- 
wardly be extended to real times on the Keldysh contour. 
For example, for the one-particle Green function one ob- 
tains 

Gi(e';0 = GJ(C';0 







-1" 


71 -C 


Go' 





where 



G, 



Go 



with the proper one-particle self-energy E. This implies 
the relation S = ^71. We have a matrix structure not 
only with respect to ^ and ^' (as in equilibrium), but also 
with respect to the Keldysh indices. How does this ad- 
ditional structure manifests itself in the Green functions 
(|13p ? For m = 1 we have two derivatives with respect to 
spinor fields, i.e. the structure of a tensor product, which 
can be made explicit by using a tensor product notation 



Gl (C',0 = (C*) 



?7— 77— 



leading to the matrix 



(18) 



Up to now we extended standard text book 
manipulations^ from imaginary time to the real-time 
Keldysh contour. We emphasize that no translational 
invariance in time is assumed. To derive the fRG flow 
equations in nonequilibrium we can follow the steps of 
Ref. \2d- In Eqs- pT|) and we replace the noninter- 
acting propagator by a propagator Gq depending on a 
parameter A G [Ap, 0] and require 



= 



A=0 



= Go 



(19) 



i.e. at the starting point A = Aq no degrees of freedom 
are "turned on" while at A = the A free problem is 
recovered. In models with infrared divergences A can 
be used to regularize the problem. In equilibrium this 
is often be achieved by implementing A as an infrared 
cutoff in momentum or energy. One of the advantages 
of the fRG approach over other RG schemes is that one 
is not restricted to these choices and other ways of in- 
troducing the parameter A have turned out to be useful 
for equilibrium problems, "^^'"^^ All that is required to de- 
rive the fundamental flow equations (|23p and (|B1[) (see 
below) are the conditions Eq. In our application 

of the nonequilibrium fRG to the steady state transport 
through an interacting quantum dot it is natural to im- 
plement A as an energy cutoff. However, such a choice 
must not be the natural one in cases where one is inter- 
ested in studying transient behavior. In this situation 
the propagator and the vertex functions in general de- 
pend on several spatial and time variables and there is 
no obvious momentum or energy cutoff scheme. Within 
the fRG several ways of introducing A can be worked out, 
compared and the one best suited for the problem under 
investigation can be identified. 

Through G^ the quantities defined in Eqs. HI]) to (fT7|) 
acquire a A-dependence. Taking the derivative with re- 
spect to A one now derives a functional differential equa- 
tion for r^. From this, by expanding in powers of the 
external sources, an infinite hierarchy of coupled differ- 
ential equations for the 7^ is obtained. Although the 
steps in the derivation are formally equivalent to Ref. [13, 
because of the real-time formulation additional factors i 
and signs appear in several places. We thus believe that 
it is helpful to present the details of the derivation for 
the present approach, which is done in Appendix A. 

In particular, for the flow of the self-energy one finds 
the expression 



FIG. 2: Diagrammatic form of the flow equation for 71 . The 
slashed line stands for the single scale propagator S^. 
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FIG. 3: Diagrammatic form of the flow equation for 72 . The 
slashed line stands for the single scale propagator S^, the 
unslashed line for G^. 



Tr 



5V(e',-;er) 



(20) 



which can be visualized by the diagram in Fig. [21 The 
trace in Eq. (j20p is meant to run over all quantum num- 
bers and time respectively frequency. In Eq. ((20)) appears 
the so-called single scale propagator (the slashed line in 
Fig.O 



= — 
dA 



(21) 
(22) 



and the quantity 72^(C'i S ') denotes the matrix ob- 
tained by keeping the indices ^ and ^' fixed. 

We thus arrive at an expression that is formally iden- 
tical to Eq. (19) in Ref. [2J. The difference appears in 
the matrix structure, which now also contains the index 
components for the branches of the Keldysh contour. To 
make this explicit, we write out Eq. (|20p with respect to 
the Keldysh structure 



a/3, A 



Tr^5^''^Sr'"^'^(r,-;e,-), (23) 



where a, (3, and v denote Keldysh indices and take the 
values =F. Apparently, the derivative of 7;!^ is determined 
by and the two-particle vertex 7^. Thus an equation 
for 7^ is required. The structure of this equation turns 
out to be similar to Eq. (21) of Ref. '2!. Here we only 
show the diagrams representing it in Fig. [3l the full ex- 
pression is given in Eq. (|B1[) . The differential equation 
for 7^ does not only contain 7;!^ - implicitly via the prop- 
agators - and 7^, but also the three-particle vertex 7;^. 
The flow of the three-particle vertex depends on the four- 
particle vertex etc. It is generically impossible to solve 
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the full set of infinitely many coupled differential equa- 
tions. In applications one has to truncate it, and this is 
usually done at order m = 2, i.e. one replaces all vertices 
with m > 2 by their initial values, which for problems 
with a two-particle interaction means — — for 
771 > 2. Even within this truncated system the remain- 
ing set of differential equations must typically further be 
approximated to allow for an analytical or a numerical 
solution,^i^iii At the end of the fRG flow, that is for 
A = 0, the 7^"° = 7m present approximations to the 
many-body vertex functions from which observables can 
be computed (for examples see Sec. IIIII) . 

We again note that a similar set of flow equations was 
derived by Jakobs and Schoeller using a diagrammatic 
approach i^^i^ 



III. APPLICATION TO QUANTUM DOTS 

The derivation of the flow equations in the previous 
section and the appendices was completely general and 
they can be used to study time dependent or station- 
ary properties of systems at finite temperature T or 
T = 0. A simple but nonetheless nontrivial application 
is nonequilibrium, stationary transport through an inter- 
acting single-level quantum dot at T = 0. For a station- 
ary situation, time translational invariance holds and one 
can Fourier transform to frequency space. Moreover, be- 
cause the quantum dot is a zero-dimensional structure, 
no additional degrees of freedom except spin have to be 
taken into account, thus considerably reducing the com- 
plexity of the problem. 

Nonequilibrium theory of single-level quantum dots 
described by the SIAM (see below) has been a major sub- 
ject of research over the past years. Several techniques 
have been developed respectively applied, for example 
statistical^-ia^, perturbative^^'S^'^^iS^^^ 
RG basedi^ii^ii^ and numerical procedures.^^ In partic- 
ular the RG based calculations were up to now restricted 
to situations away from the Kondo regime, either by 
looking at features in the mixed-valence statei^ or in 
strong magnetic fieldsJ^ii^ Despite this tremendous 
ongoing effort, relatively little is known about nonequi- 
librium properties, even in the stationary state, of the 
model. On the other hand, the very detailed knowledge 
of the equilibrium properties^ makes it possible to 
interpret certain features or even dismiss approximations 
on the basis of this fundamental understanding. 



A. Single impurity Anderson model 

The standard model to describe transport 
through interacting quantum dots is the SIAM.^^ 
Experimentally;^'^ the dot region is attached to two 
external leads, and the current is driven by an applied 
bias voltage Vb ■ Furthermore, the filling of the dot can 
be controlled by a gate voltage Vg ■ To keep the notation 



simple we use units with e 
reads 



1. The Hamiltonian 



H 



ka-l kal kal 

kal 



E 



+ 



kal 



VadU^ + U 



H.c. 



(24) 



in standard second quantized notation {n„ = d\,d^). As 
usual, k is the wave vector of the band states in the leads 
and a the spin. In addition, the index I = L, R dis- 
tinguishes the left and right reservoirs which can have 
different chemical potentials /i/ through an applied bias 
voltage Vb — /^l — fJ-B- Since we do not include a mag- 
netic field lifting the spin-degeneracy of the spin up and 
down level, the one-particle Green function and the self- 
energy are spin independent and we suppress the spin 
index. We assume the dispersions e^r^^ and hybridiza- 
tions Vj^^i between dot and the left and right leads to be 

identical, and Vf: = V/\/2 to be fc- independent. Equation 
p4|) is written such that for Vb = we have particle-hole 
symmetry. 

The interaction in the Hamiltonian is reduced to the 
dot site only. We can thus reduce the problem to a zero- 
dimensional one by "integrating out" the leads. For the 
flow equations we need the propagator Gd.o at C/ = 0. 
To this end we use the Dyson equation to obtainii^i^ 



where 

Gd,o 



Lu - Vg 






-uj + Vg 



- S 



lead 



(25) 




V 



-'lead — 



— T 

2N ^ 



kl 



d,0 ^ 




++ ; 




^kl' 


^kl 


G+- 

kl 


G+ 

kl 



(26) 



To further evaluate Eq. (|26|) , we have to insert the Green 
functions of the free electron gas, given by 

1 



^ki^-) 



) - ej: + y.1 + iS 
-2iTrf{ef:)S{uj - ej:- 



^ki^-) 



2iT:f{e^)5{uj- + III) , 
-2i7r/(-£^)(5(w -e^ + ^n) , 

into Eq. This leads to 

LO-VG-iT [1 - f^iu) - fniLu)] 



{uj-VGr + r^ 



(27) 

(28) 
(29) 
(30) 

(31) 
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Gdii^) = -[Gjfli^)]* , (32) the leads and T = tt\V\'^Nf, with Np the density of states 

r [/^(ll!) + fji(uj)] Fermi level of the semi-infinite leads, represents the 

G^Q (oj) — i—^ y -|2'_|_ p2 ' ^^"^^ tunnel barrier between the leads and the impurity. 

_i ,r cj) + tj)] Finally, the quantities we want to calculate are the cur- 

d,Q\^) — * (w — Vg)'^ + ' rent J through the dot and the differential conductance 

G — dJ/dVB- For the model ([M]) the current is given 

where fi{±Lo) — f — /i;)) are the Fermi functions of by^'^^ 



J - ^{Jl + Jn)^§^Jde [hie] ~ fnie)] [G+~{e) - G^+(e)] , (35) 

where we already used that the left and right couplings are identical, i.e. Tl = Tji = r/2, and summed over both 
spin directions. The interacting one-particle Green function of the dot is denoted by Gd and Jl/r are the currents 
across the left and right dot-lead contacts respectively. Equation psp is written in a somewhat unusual form, not 
employing the relation G'^~ {e) — G'^~^{e) = G^(e) — G"^{e) = —2iTipd{e), where pd denotes the dot's one-particle 
spectral function. The reason for this will be explained below. 

Another quantity of interest is A J — Jl — Jr^ Obviously, since no charge is produced on the quantum dot, 
A J = in the exact solution. Using again the results of Refs. andlH, the expression for A J becomes 

TT 7 A(CJ) ^ ^ 

— oo 

A(w) = \uj -VG + iT[l- F{uj)]-Y.— {uj)\^ +\rF{-Lj)+^rrii:+-{Lj)\ [rF(w) - !3ml]"+(a;)] 

F{lo) = fL{uj) + fR{uj) (37) 



where we used that S (oj) and [uj) are purely imag- 
inary. Depending on the type of approximation used 
A J — might either hold for all parameters^^'^^ or not.^^ 
We note that fulfilling A J = is, however, not suf- 
ficient for an approximation to provide reliable results. 
E.g. the self-consistent Hartree-Fock approximation ful- 
fills A J = 0, but nonetheless does not capture the correct 
physics even in equilibrium .~ 



B. Lowest order approximation 

Before studying the coupled system for E'^ and 7^ we 
begin with the simpler case where we replace 7^ on the 
right hand side of Eq. ([23]) by the antisymmetrized bare 
interaction and consider only the flow of S. The only 
nonzero components of the bare two-particle vertex are 
those with all four Keldysh indices being the same (a = 



d_ 

dK' 



Within this approximation the self-energy is always time 
or frequency independent, and no terms off-diagonal in 
the Keldysh contour indices are generated. It leads to at 
least qualitatively good results in equilibrium,— and has 
the additional advantage that the flow equations can be 
solved analytically. 

As the last step we specify how the parameter A is 
introduced. Since we are interested in a stationary situ- 
ation, i.e. the propagators only depend on the time dif- 
ference t — t' , all equations can be transformed into fre- 
quency space and one natural choice is a frequency cutoff 
of the form 



GUco)^ei\iu\-A)GdA^) 



(40) 



with Aq — > ooi^ Evaluating S by means of the Morris 
lemmaSi"^ results in 



(38) 



S^{lj)^S{\oj\-A) 



1 



Gd,o{^) -S^ 



(41) 



With this replacement Eq. (|23p reduces to 

-E^T.A ^ ±iu f ^5TtA(^) . (39) 
J 27r 



A straightforward calculation permits us thus to rewrite 
Eq. ([391 as 
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^IM _ v±±, 

A(^) ^ 



(42) 



where 



1 



Finally, the initial condition for the self-energy is 
hm t^'> = 0^ 

Aq — ^oo 



1. Equilibrium 

We now focus on T = 0. In a first step we discuss the 
equilibrium situation, that is Vb — I^-l ~ I^r ^ 0. Then 
we obtain the decoupled system 



dA TT [(A ± ir)2 - (Vb ± S=FT, 



(43) 



which can be solved analytically. We first note that with 
1]++'^ = — [E~~''^] both equations are equivalent. For 
Y,^^'^ we obtain with the definition = Vq + T,^^'^ 
the solution 



(44) 



where Jn and F„ are the Bessel functions of first and sec- 
ond kind. The desired solution of the cutoff free problem 
is obtained by setting A = 0, i.e. 



(45) 



which is precisely the result Eq. (4) obtained by Ander- 
gassen et ai— It is, however, important to note that in 
this work, based on the imaginary-time formulation of 
the fRG, the differential equation has a different struc- 
ture. It is real and has a positive definite denomina- 
tor. Thus, while the solutions at A = are identical 
for the imaginary-time and real-time formulations, the 
fiow towards A = will show differences. As we will see 
next, the complex nature of the differential equation (j43p 
can lead to problems connected to its analytical struc- 
ture when attempting a numerical solution. For small 
U/T no particular problems arise. As an example the 
result for the flow of S as function of A for U/T = 1 
and Vg/F = 0.5 obtained with a standard Runge-Kutta 
solver is shown in Fig. U) Consistent with the analytical 




-0.2 



10^ 10° 10^ 10^ 10^ 

A/r 



FIG. 4: Flow of E~''^/r with A/F for U/T = 1, Vg/F = 0.5, 
and Vb = 0. The full curve shows the real part, the dashed 
the imaginary part of E 




-15 



-20 



— A=oo^A=o ■ 
■-■ A=0->A=oo 
^ X exact solution- 



10^ 10'^ 10'^ 10° 
A/F 



10^ 10^ 10^ 



FIG. 5: (color online) Flow of S—'^^/F with A/F for U/T = 
15, Vg/F = 6, and Vb = 0. The full and dashed curves show 
real and imaginary part obtained from the integration A — 
CO — > A = 0, the dashed-dotted and dotted curves real and 
imaginary part obtained from an integration A = — > A = oo, 
using the solution from (|45|l as initial value for E The 
crosses denote the analytical solution (|44)) . 



solution Eq. (j45|) . the imaginary part (dashed line) goes 
to zero as A ^ 0, while the real part (solid line) rapidly 
approaches the value given by formula (j45p . 

However, for larger values of U/T the numerical so- 
lution becomes unstable in a certain regime of Vq. A 
typical result in such a situation is shown in Fig. [5] The 
different curves were obtained as follows: The full and 
dashed ones from the numerical solution starting with 
E = at Aq — > oo, the dash-dotted and dotted 
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FIG. 6: (color online) Imaginary part of Ap determined from 
(|44|) and Va + S '^^ = Ap + ir for different values of U and 
Vb = as function of Va- For U /V = 15 and 20 there exist 
an interval [Vg, Vg] where SJmAp > 0, while for small U or 
Vg [V§, Vq] we always have SmAp < 0. 



by integrating the differential equation (|43p backwards 
from A = with the correct solution for A = as given 
by formula (|45p as initial value. The crosses finally are 
the results from the analytical solution Eq. (|44|) . Ev- 
idently, there exists a crossing of different branches of 
solutions to the differential equation for A/F w 1 and 
the numerical solution with starting point A = oo picks 
the wrong one as A ^ 0. The reason for this behavior 
is that for large U there exists a certain Vq such that 
Vq + E^^'^^p = Ap + iF with real Ap, resulting in a pole 
in the differential equation (|43p . For Vg Vq this pole 
does not appear for real Ap, but as shown in Fig. [6] 3m Ap 
changes sign at Vq, which in turn induces a sign change 
on the right hand side of the differential equation, lead- 
ing to the behavior observed in Fig. [51 There also exists 
a second critical value Vq such that for Vg > Vq we find 
SmAp < and the instability has vanished again. 

Obviously, this instability limits the applicability of the 
present approximation to sufficiently small values of U. 
This is different from the imaginary-time approach by 
Andergassen et al.^ where this simple approximation 
leads to qualitative correct results even for values of U 
significantly larger than F. 



2. Nonequilibrium 

We now turn to the case of finite bias voltage Vb- As a 
typical example, the flow of S~~'^ for U/T = 1 (full and 
dashed curves) and 5 (dashed-dotted and dotted curves) 
for Vg/T = 0.5 at Vb/T^O (equihbrium) and Vb/T = 1 
is shown in Fig.[71 Since the results for 1]++'^ are related 
to those for S~'^ by T,++'^ = - [T.—^^]* we do not 
show them here. The Vb dependence of the curves for 
SRe [E (thick lines) looks sensible. For Vb ^ an 





• h ■ - ' 


' / 


- 










1 — u/r=i 


Vg/r=o _ 


/ _ . u/r=i 


VB/r=i . 


.- u/r=5 


VB/r=o . 


... u/r=5 

1 1 1 1 


VB/r=i . 
1 



I c;! I I I I I 

-? -1 n 1 ? T 

10 10 10 10 10 10 10 

A/r 



FIG. 7: (color online) Flow of T,~~-^/T with A for U/V = 1 
and 5 for Vg/F = 0.5 at Vs/F = and 1 (thick curves: 
real part; thin curves: imaginary part). The curve for 
SimE— •'^/r at U/T = 1 and Vs/F = 1 (thin dashed line) 
lies on top of the corresponding zero bias curve and is thus 
not visible. 

imaginary part of order is generated in the flow which 
does not vanish for A ^ (see the thin dotted line). 
Causality requires that the relation 

J:"iuj) + E++{u;) = - [E-+(c^) + E+-(cj)] (46) 

must hold for the exact solution. Because of S ''(w) = 
E"' {uj) = 0, the finite imaginary part of S"" leads to a 
breaking of the condition (|46)) to order at the end of 
the fRG flow. This is consistent with the fact that by ne- 
glecting the flow of the vertex terms of order are only 
partially kept in the present fRG truncation scheme. To 
avoid any confusion we emphasize that our RG method 
is different from any low-order perturbation theory. The 
weak breaking of causality can also be understood as a 
consequence of our approximation leading to a complex, 
energy-independent self-energy: The off-diagonal com- 
ponents, being related to the distribution functions for 
electrons and holes, respectively, in general have differ- 
ent support on the energy axis. The energy independence 
makes it impossible to respect this structure here. 

For our further discussion the order [/^ violation of 
Eq. means that we may not rely on relations like 
Eqs. (Ell-dTni) but have to work with G"'', thus the 
somewhat unusual formula (|35[) . A naive application of 
E^ = E~- - E-+ and use of Gj' - G~+ = 2i3mGf 
would have led to unphysical results. That working with 
G'^^ — G^^ is still sensible can be seen from a straight- 
forward evaluation leading to 

G+-{u)-G-^+{u)^ (47) 



|w - Vg + «r [1 - F{oj)] - E— 1^ + F2F(w)F(-w) 

F{u;)= fL{Lo) + fni^) 

which is purely imaginary with a definite sign. Insert- 
ing the expression (|T7)) into the formula (pS)) . one can 
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FIG. 8: (color online) Conductance normalized to Go ~ 2e^ /h 
as function of Va for U/T = 2 and several values of the bias 
voltage Vb- 

calculate the current and thus the conductance. Since 
we are at T = 0, an explicit expression for the current 
of the cutoff free problem (at A = 0) can be obtained 
by noting that with /i^ = Vb/2, ^.r — ~VbI'2' one has 
- = e(VB/2 - Icc-I) and F{±lo) = 1 for 

Lo e [-Vb/2, Vb/2], which leads to 

2 ^''^^ 

J — / duj- TT 

-Vb/2 

with the abbreviations 

Vg = Vc + ^eE— , (49) 

r* = + (SJmE— )2 . (50) 

Equation pS)) for the current is equivalent to the non- 
interacting expression but with renormalized parameters 
Vq and r*, which depend on the interaction as well as 
the bias and gate voltage. 

An example for the differential conductance as func- 
tion of Vg obtained from Eq. for U/T = 2 and sev- 
eral values of Vb is shown in Fig. [51 where Gq — 2e^//i 
(after reintroducing e and h). Increasing Vb leads, as 
expected, first to a decrease of the conductance close to 
Vg — and later to a splitting of order Vb ■ Since we will 
discuss a more refined scheme including parts of the flow 
of the two-particle vertex next, we do not intend to dwell 
too much on the results of this simplest approximation. 
We note in passing that for Vg = due to particlc-holc 
symmetry we obtain from the differential equation (|42p 
that S^^ = independent of U. Consequently the cur- 
rent J calculated via Eq. (|48)) and the conductance are 
independent of U, too, and given by the corresponding 
expressions for the noninteracting system. As we will see 
in the next section, this deficiency will be cured by the 



FIG. 9: (color online) Flow of E'-'^/F with A for U/T = 1 
(left panel) and U/T = 15 (right panel) at Vg = U/2 and 
Vb = 0. The full curves show the real part, the dashed the 
imaginary part of E~~'"^. The stars at the vertical axis denote 
the values as obtained from the imaginary-time fRG.— 

approximate inclusion of the vertex flow. In the present 
approximation the current conservation A J — holds for 
all parameters as E-+ = S+- = [cf. Eq. 



C. Flowing vertex 

A more refined approximation is obtained when we in- 
sert the flowing two-particle vertex 7^ as given by ex- 
pression (|B2[) in the calculation of the self-energy Eq. 

By this we introduce an energy-dependence of the 
self-energy<^ However, because the size of the resulting 
system of differential equations becomes extremely large 
if the full frequency dependence is kept (for a discussion 
on this in equilibrium see Ref . Ilil ) , we only keep the flow 
of the frequency independent part of the vertex, an addi- 
tional approximation which has successfully been used in 
equilibriumi^S As a consequence we again end up with a 
frequency independent S^. The resulting expression for 
the self-energy (see Appendix B) is 

uj=±A 

(51) 

where U"'^^'^'^ is the flowing interaction given by the 
expression (|B3p . As has been observed by Karrasch et 
al.,^'^ this approximation leads to a surprisingly accurate 
description of the transport properties in equilibrium. In 
particular it is superior to the lowest order approximation 
including only the bare vertex. 



1. Equilibrium 

We again begin with the discussion of the solution to 
Eq. (|?T|) in equilibrium. Results for the flow of S are 
presented in Fig.Ofor U/T = 1 (left panel) and U/T = 15 
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FIG. 10: (color online) Flow of S" '^/F and E'+'^/F with 
A for U/V = 1 (left panel) and U /V = 15 (right panel), 
Vg/F = U/2 and Vb/F = 1. The full curves show the real 
part, the dashed the imaginary part of E , the dot-dashed 
the imaginary part of E""'"''^. The real part for the latter is 
zero. 

(right panel) for Vg = U/2. Since = - [S^--^]* 

only one component is shown. The stars in Fig. [9] denote 
the solutions of the imaginary-time equations taken from 
Ref.m. Note that for U/T = 15 and Vfe/F > 6 the simple 
approximation Eq. (|42p showed an instability, while with 
the flowing vertex the system is stable even for these 
large values of U and reproduces the correct equilibrium 
solutions for A Oi^ The reason for this is that the flow 
of the vertex reduces the resulting effective interaction 
below the critical value in the instability region.— 

2. Nonequilibrium 

For the same parameters as in Fig. [S] we present the 
resulting flow with finite bias Vs/r = 1 in Fig- UHl In 
addition to the curves for real (solid lines) and imaginary 



part (dashed lines) of S a third curve is displayed, the 
imaginary part of E (dashed-dotted lines), which now 
is generated during the flow. Note that SfteS ^'^ = 
and S'' = [E '"'^] . Furthermore, we always find 
3mS-+''^ < 0. For U/T = 15 (right panel in Fig.[TOl) we 
have rescaled ^toE by a factor 10^ to make it visible 
on the scale of S Since S"'''^ is a complex energy 

independent quantity Eq. is again not fulfilled. We 
note that the error is still of order , but for fixed Vg 
and Vb it is significantly smaller than in the simplest 
truncation scheme discussed above. 

The energy-independence of the self-energy allows to 
derive an analytical expression for the current at T = 
similar to Eq. (|i5)) . which due to the appearance of 
now becomes 

J=Iy^^^' arctan (^^^^^ j (52) 
with as in Eq. ^ and 

f = F - 3ml]"+ > F 

F* = ^f2 + (QmE— )2 ; 

where T,"^ is taken at A = 0. Thus, the only change to 
the expression is a formal replacement F F in J/F. 
Equation (|52[) is of the same structure as for the nonin- 
teracting case with Vg and F replaced by renormalized 
parameters. However, the two self-energy contributions 
S and E ^ enter distinctively dilferent in the expres- 
sion for the current. While Q^toS solely plays the role 
of an additional life-time broadening, StoE ^ directly 
modifies the tunneling rate both in the prefactor of J 
and in the expression for the life-time broadening. 

A problem occurs when using the results of the present 
approximation in Eq. (|36p. leading to 



F f , 3mE-+ [1 - F(uj)] , , 

AJ = 2- / du: 5 ^ "—^ . 53 

71" 7 |w - Vg + ?T [1 - F(cj)] - E— r + [FF(-cj) + 3mE+-] [FF(cj) - 3mE-+] 



The requirement AJ = is only fulfilled for Vg = 0, 
because then E =0 and the integrand is asymmetric 
with respect to uj. Thus our approximation of an energy- 
independent flowing vertex violates current conservation 
for Vg 7^ in nonequilibrium. We verified that A J ~ U'^ 
which is consistent with the fact that not all terms of 
order C/^ are kept in our truncated fRG procedure. How 
does A J behave in the limit Vb 0? To see this we note 
that, because 3mE does not depend on the sign of Vb 
and furthermore goes to zero as Vb 0, 3mE ^ 



Vg. Consequently, A J ^ and hence the violation 
of current conservation vanishes in the linear response 
regime Vg ^ 0. 

In Fig. [11] we show the current at Vg = as function of 
Vb for t//F = 1, 6 and 15. With increasing U the current 
for intermediate Vb is strongly suppressed. In addition 
there occurs a structure at low Vb , which turns into a re- 
gion of negative differential conductance with increasing 
\J . The appearance of such a shoulder in the current was 
observed in other calculations as welli^'^'^ However, 
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FIG. 11: (color online) Current normalized to Jo = Go 7 
(after reintroducing e and h) as function of Vb for U/T = 1, 
6 and 15 and Vg = 0. For U /T — 15 we find a region of 
negative differential conductance in the region IVb/FI ~ 0.5 
(c.f. Fig.[T2j. 



FIG. 13: (color online) Differential conductance G as function 
of Vg for different values of Vb and U/F = 1 (upper panel) 
and U/r = 15 (lower panel). Note the extended plateau at 
Vb ~ for U/T — 15, which is a manifestation of the pinning 
of spectral weight at the Fermi level. 
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FIG. 12: (color online) Differential conductance G as function 
of Vb for Vg ~ and various values of U. For [//F > 5 a 
distinct minimum around Vb /F ~ 0.5 appears. 



whether the negative differential conductance we find for 
still larger values of U (c.f. Fig. fT^ is a true feature of the 
model or rather an artifact of the approximations used is 
presently not clear and should be clarified in further in- 
vestigations. However, negative differential conductance 
has also been observed in a slave-boson treatment of the 
modelj^ 

Keeping Vg = fixed, we can calculate the conduc- 
tance G — dJ/dYs as function of Vb for different values 
of U. The results are collected in Fig. [T^l In contrast 
to the simple approximation without fiow of the vertex, 
the conductance is now strongly dependent on U, except 
for Vb — 0, where due to the unitary limit at T = we 
always find G — Gq- As already anticipated from the cur- 
rent in Fig. 1111 a minimum in G starts to form around 
Ve/r « 0.5 for U/T > 5, which is accompanied by a 
peak at Ve/r w 2. A similar behavior in the conduc- 
tance was observed in a perturbative treatment, which 
in contrast to our current approximation involves the full 



energy-dependence in the self-energy. This at least quali- 
tative agreement - we of course cannot resolve structures 
like the Hubbard bands with an energy independent self- 
energy - again supports our claim that despite the viola- 
tion of the relation (|46|) we can obtain reasonable results 
from G"^. 

We finally discuss the variation of the conductance 
with Vg for fixed U and Vb- We again emphasize, that 
for Vg 0, AJ = only holds to leading order in U. 
In Fig. [13] we present the curves for two different val- 
ues of U, namely U/T = 1 (upper panel in Fig. [75]) and 
U/r — 15 (lower panel in Fig. [T3)) . In the former case, 
the variation of G with Vb is rather smooth, as is to 
be expected from the current in Fig.[TlJ For large U, we 
observe an extended plateau at zero bias, which is a man- 
ifestation of the fact that in the strong coupling regime 
a pinning of spectral weight at the Fermi energy occurs. 
This feature is also observed in the imaginary-time fRG 
as well as in NRG calculations."^^ Increasing Vb quickly 
leads to a similarly extended region of negative differen- 
tial conductance, which, assuming that this result is a 
true feature of the model, therefore seems to be linked 
to the "Kondo" pinning. We note that it is unlikely that 
the appearance of the negative differential conductance 
is related to the breaking of current conservation at order 
as it also appears for Vq = where A J = 0. For large 
Vb multiple structures appear in G, which are related to 
the energy scales Vb and U. 



IV. SUMMARY AND OUTLOOK 

Starting from a generating functional proposed by 
Kamenee\^ we have derived an infinite hierarchy of dif- 
ferential equations for the vertex functions of an interact- 
ing quantum mechanical many-body system in nonequi- 
librium (see also Refs. [25l and[26l). The major difference 
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to the imaginary-time formulation comes from the use of 
Keldysh Green functions. The indices + and — referring 
to the upper and lower branch of the (real-)time contour 
lead to an additional matrix structure to the problem. 
Our formulation is sufhciently general that it allows to 
treat bosonic and fcrmionic models with or without ex- 
plicit time dependence and at T > 0. Since the fRG 
leads to an infinite hierarchy of coupled differential equa- 
tions, one has to introduce approximations, at least a 
truncation at a certain level. Typically one neglects the 
flow of the three-particle vertex. As has been demon- 
strated in Ref. [2J for the imaginary-time fRG one can 
solve the remaining system of flow equations for simple 
models like the SIAM numerically, thus keeping the full 
energy-dependence. Due to the fact that the vertex func- 
tion carries three continuous frequency arguments in ad- 
dition to the discrete quantum numbers of the system, 
such a calculation can become computationally quite ex- 
pensive. 

To reduce the numerical effort, further additional ap- 
proximations can be introduced. A particularly impor- 
tant and successful one is obtained by neglecting the en- 
ergy dependence of the vertex functionSf^^i^S which al- 
ready leads to a surprisingly accurate description of local 
and transport properties of interacting quantum dots in 
the linear response regime. 

We applied the nonequilibrium fRG to the SIAM with 
finite bias voltage in the stationary state. It turned out, 



that for the simplest approximation where only the flow 
of the self-energy is kept, the analytic structure of the 
differential equation leads to problems in the numerical 
solution. In addition, this approximation leads to a vio- 
lation of the causality relation (|46p to order . The first 
problem was resolved by including the two-particle vertex 
in the flow at least up to the largest interaction consid- 
ered here {U/T = 15). At the present stage this was for 
computational reasons done by assuming it to be energy- 
independent, yielding again an energy-independent self- 
energy. Although this approximation also violates Eq. 
PS)) to order for a fixed Vg and Vb the error is sig- 
nificantly smaller compared to the simplest scheme. We 
were able to obtain reasonable expressions and numer- 
ical results for the current and the conductance using 
the functions G°'^{uj) instead of G^{lo) in the current 
formula. We reproduced nonequilibrium features of the 
current and differential conductance known from the ap- 
plication of other approximate methods to the SIAM. 

In the more advanced truncation scheme and for Vq 
0, the current conservation A J = only holds to lead- 
ing order in U . This defect can be traced back to the 
energy independence of the two-particle vertex, leading 
to finite, but energy-independent S and S"' . By a 
close inspection of formula (|36p . however, one can see 
that this deficiency can for example be cured by assum- 
ing a coarse-grained energy dependence of the form 



J 



— < < — : F{^) = 1 



^ < w : F{uj) = 



^ QmT,+-^^ = , !3mS-+'^=const. 



^ 3mS+"-^ = -3mi;"+'^ = const. 



Q?mi;"+-^ = , 5mE+-'^==const. 



(54) 
(55) 
(56) 



and corresponding energy dependencies for E and 
E++. Here the const, might depend on A but not on oj. 
Such an approximation will quite likely also reestablish 
causality Eq. (|46p . Work along this line is in progress. 

Including finite temperature, magnetic field etc. in the 
calculations of transport is easily possible, as well as the 
extension to more complicated "impurity" structures (see 
e.g. Ref. [39I) . The latter aspect is typically a rather cum- 
bersome step for other methods like perturbation theories 
or in particular numerical techniques. 

Our present stage of work should mainly be seen as a 
"proof of principle" . Already for the SIAM there remain 
a variety of fundamental things to do; like the implemen- 
tation of the full energy-dependence in the calculations, 
which then should cure the violation of charge conser- 
vation and causality; or implementing a reasonable "cut- 



off scheme" for time-dependent problems off equilibrium, 
e.g. to calculate current transients. Despite its current 
limitations we expect the real-time nonequilibrium for- 
mulation of the fRG to be a useful tool to understand 
nonequilibrium features of mesoscopic systems, just like 
the imaginary-time version in the calculation of equilib- 
rium properties. 
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APPENDIX A: DERIVATION OF THE FLOW EQUATIONS 

Although the following derivation is mainly equivalent to the one presented in Ref. [23, the appearance of the factor 
i in the real-time formulation of the generating functional leads to some changes in signs and prefactors i in the 
equations. It thus appears to be helpful for the reader to repeat the derivation. 

As a first step we differentiate W^'^ with respect to A, which after straightforward algebra leads to 



dK 



c,A 



C Tr (q^G°'^) +<Tr (q 



SfjSrj 



(Al) 



Considering and </> as the fundamental variables we obtain from Eq. (|16p 



m, m = -^w^' ({r }, m) - ( 4> 



dK 



dK 



Applying the chain rule and using Eq. (|A1[) this leads to 



where the last term in Eq. (|16p cancels a corresponding contribution arising in (jAip thus a posterior justifying the 
inclusion of this term. 

Extending the well known relationii between the second functional derivatives of T and to nonequilibrium we 
obtain the functional differential equation 



dK 



-C Tr ( Q^G"'^ 



Tr(Q^V|i(^^G"'^• 



(A2) 



where Vj-'^ stands for the upper left block of the matrix 



V 



-c 



G' 



O.A 



GO- 



(A3) 



and the upper index t denotes the transposed matrix. To obtain differential equations for the 7^ which include 
self-energy corrections we express Vj, j, in terms of instead of G"-^. This is achieved by defining 



7/ ^ A 



and using 



G' 



0,A 



C7f 



(A4) 



which leads to 



^ _^ Tr (q^G°-^) -KCTr [G^Q^V|'^(r^, G 



(A5) 



with 
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It is important to note that ^ as well as j^g^ and j^g^ are at least quadratic in the external sources. The initial 
condition for the exact functional differential equation (|A5[) can either be obtained by lengthy but straightforward 
algebra, which we are not going to present here or by the following simple argument: at A = Aq, G^'^" = (no 
degrees of freedom are "turned on" ) and in a perturbative expansion of the 7,'^^° the only term which does not vanish 
is the bare two-particle vertex. We thus find 



(A7) 



Expanding m Eq. (jA6p in a geometric series 



G' 



* 5^T 



(A8) 



and r in a Taylor series with respect to the external sources 



(A9) 



one can derive an exact infinite hierarchy of fiow equations for the 7^. For example, the flow equation for the single- 
particle vertex 71 (the self-energy) can be obtained by taking the expansion (jA8l) up to first order in U^^^^ inserting 
it into (|A5p . replacing on both sides by the expansion (jA9p and comparing the terms quadratic in the external 
fields. This procedure leads to the expression ((20|l . 



APPENDIX B: FLOW EQUATION FOR THE TWO-PARTICLE VERTEX 



To obtain the flow equation for E'^ we have expanded Eq. \k&\ in a geometric series up to first order and F'^ (|A9|) 
in the external sources up to m = 2. To find the flow equation for the vertex function, 7^, we have to proceed one 
step further, that is expand l|A6|) up to the second order and F'^ up to m = 3. After comparison of the terms with 
the same power in the fields we obtain an expression similar to the equilibrium one (see Eq. (15) of Ref. 

+ 7"''^''^^ (fi , ^3 ; a , UW^-''^'' i^'i , ; 6 , 6 ) 
- 7''^^'''^(e2, ^4; ^ui3h''"''^\e3,ei;U, 6) 



where the contribution involving 7;^ was replaced by its 
initial value, i.e. 7;^ — — 0. The main difference 
to equilibrium is the presence of mixed contractions for 
the Kcldysh indices of the matrices G^ and and the 
four-index tensors 7^. 

Up to this point the derivation of the fiow equation 
for the two-particle vertex was as general as the discus- 
sion in Sec. |lll We now consider the application to the 
steady state current through a quantum dot described by 
the SIAM presented in Sec. IIIII We thus go over to fre- 
quency space and introduce the cutoff A as in Eq. (|40)) . 
The single-particle quantum number in ^ is given by the 
spin index. We can take advantage of spin and energy 
conservation which implies that the lower indices of G 



and S are equal, that is ^3 = ^3 and ^4 = ^4. This does 
not hold for the Keldysh indices. The integrals over UJ31 
and as well as the sums over a'^ and <t'^ can then be 
performed trivially. Using the Morris lemma^^ we can 
rewrite the matrix product of G and S as 

1 



G^iu)S^i^') ^ -6{\uj'\ - K)Gi{u^)Gi{J) , 



with 



Grf,o(w)] ^ -t^{Lo) 
The (5-function can be used to perform another of the 
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frequency integrals, and the remaining one can be eval- vertex. Using the spin independence of the Green func- 
uated because of energy conservation of the two-particle tion for zero magnetic field, this leads to 



47r ^ ^ 

a;=±A 0-3 ,0-4 fj,,u,p,ri 



-GT\o.)GT\u^) 



' cr^ ,(T4;cri ,(73 ' (73 ,cr2 ;cr4 ,iT2 ' ct^ ,(73 ;cri ,(74 ' (74 ,cr2 ; crs ,cr2 



' ,(74 ,0-3 ' 173 ,(7^ ;(74,C72 



J. 



p.a.ri5^A 



2,i73;f7i,Cr4 '(74,(7-^ ;(73, (72 



(B2) 



Comparing (jB2[) with Eq. (20) in Ref. [39l . were similar approximations were made in equilibrium, we see that we have 
two more terms because of the Keldysh indices. In the first term appears C^' , while its transpose G^^/' enters 
everywhere else. 

Using the antisymmetry of 7^ in the spin indices and spin conservation we can further simplify Eq. (|B2p by 
introducing the flowing interaction U"^^^'^ defined as 

Q/37(5,A — X , X , Tjap-tS.K _ X , X , TT0afS,A 
I(7[,cr'^;ai,a2 ~ "ctJ ,cri "cr^ ,(T2 "cr^ ,cri "o-;^ ,0-2 

This leads to the flow equation 

J^jjaff-fS.A _ ^ ^ f Qprl.A^_^•^Ql^^J.,A^^^^ ^-^a(3piy,A^rifi'y5,A ^ ^(3ap:^,Ajjfirj-fS,A-^ 



uj=±A p,i^p,v 



-G^^^'^ (lu)G''^'^ [uj) ^2U°'^'^^''^UP^'^^'''^ — U'^f^'^^'^jjI^P'^^'^ ~ jjtJ-o--t<n,AjjppvS,A 

_^t2jj-ap^i',Ajj-pf3riS.A jj-apjiy,Ajj-f3prj5,A jjpa^v^AjjpfirjS.A 

jjp(i^ri.Ajjapv5.A jj-pf^jv.AjjaprjS^A' 



(B3) 



From the initial value of 7^ at A = Aq — * cxd given in while all other components are zero. 
Eq. (081) wc can read off as the initial value for [/"/3f'5-Ao 



U 



Oi.aaa.A 



° = aiU , 
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